Data collected from sensors are always degraded by physical blur and noise to some extent—it is impossible to build imaging instruments that produce arbitrarily sharp images uncorrupted by measurement noise. Digital image reconstruction is increasingly being used to enhance images by removing noise and blur. The role of digital image reconstruction, performed by software, hardware or a combination thereof, is to process the raw input data so that the resultant images provide the best possible approximation of the true underlying images. The reconstructed images are then presented to viewers or to machine image analysis tools in order to extract the information of interest. Image reconstruction therefore reveals information which is present, but may be hidden in the raw data by blur and noise. In image reconstruction, the main challenge is to prevent measurement errors in the input data from being amplified to unacceptable artifacts in the reconstructed image.
In many existing approaches to image reconstruction, images are deblurred and denoised by replacing the image intensities in pixels by weighted averages of intensities in neighboring pixels. Noise is reduced if the weights are all positive, because positive noise fluctuations cancel negative ones. But this also causes the signals in the pixels to be averaged, which results in added blur. This reconstruction bias can be mitigated by using a mix of positive and negative weights. Such a combination of positive and negative weights can also be used to undo the physical blur inherent in the raw data. However, weights of mixed signs can also amplify the noise, sometimes very significantly. The goal of image reconstruction, therefore, is to find the optimal weighting scheme, which maximizes noise reduction while minimizing the bias caused by signal averaging. If possible, image reconstruction should also reduce the original physical blur to the extent possible.
The bias introduced by signal averaging can be assessed by expanding the signal in a Taylor series around a pixel. To be specific, a two-dimensional example is described herein, but the formulation can be generalized trivially to any dimensionality. Let (x, y) be the position of a reference pixel that is to be denoised. The signal in the neighborhood of (x,y) can be expanded as
                              I          ⁡                      (                                          x                +                r                            ,                              y                +                s                                      )                          =                                            I              0                        +                                          ∂                I                                            ∂                x                                              ⁢                      |            0                    ⁢                      r            +                                          ∂                I                                            ∂                y                                              ⁢                      |            0                    ⁢                      s            +                                          1                2                            ⁢                                                                    ∂                    2                                    ⁢                  I                                                  ∂                                      x                    2                                                                                ⁢                      |            0                    ⁢                                    r              2                        +                                                            ∂                  2                                ⁢                I                                                              ∂                  x                                ⁢                                  ∂                  y                                                              ⁢                      |            0                    ⁢                      rs            +                                          1                2                            ⁢                                                                    ∂                    2                                    ⁢                  I                                                  ∂                                      y                    2                                                                                ⁢                      |            0                    ⁢                                    s              2                        +            …                                              (        1        )            
where the subscript 0 designates the value at the reference point (x,y), r is the displacement from x, and s is the displacement from y. The first derivatives of I are the components of the image gradient, and the second derivatives are the components of the curvature matrix. Additional terms in the Taylor expansion—not shown explicitly—include higher order polynomials in r and s with coefficients that are higher partial derivative tensors of I at (x, y).
Signal averaging is achieved by multiplying I (x+r, y+s) by a filter g (r,s) and integrating over r and s. The result is the averaged signal
                                                                                                                 J                    ⁡                                          (                                              x                        ,                        y                                            )                                                        )                                =                                ⁢                                  ∫                                                            I                      ⁡                                              (                                                                              x                            +                            r                                                    ,                                                      y                            +                            s                                                                          )                                                              ⁢                                          g                      ⁡                                              (                                                  r                          ,                          s                                                )                                                              ⁢                                          ⅆ                      r                                        ⁢                                          ⅆ                      s                                                                                                                                              =                                ⁢                                                                                                    I                        0                                            ⁢                                              μ                                                  (                          0                          )                                                                                      +                                                                  ∂                        I                                                                    ∂                        x                                                                              ⁢                                      |                    0                                    ⁢                                                            μ                      x                                              (                        1                        )                                                              ⁢                                                                  ∂                        I                                                                    ∂                        y                                                                              ⁢                                      |                    0                                    ⁢                                                            μ                      y                                              (                        1                        )                                                              +                                                                  1                        2                                            ⁢                                                                                                    ∂                            2                                                    ⁢                          I                                                                          ∂                                                      x                            2                                                                                                                                ⁢                                      |                    0                                                                                                                                          ⁢                                                                            μ                      xx                                              (                        2                        )                                                              +                                                                                            ∂                          2                                                ⁢                        I                                                                                              ∂                          x                                                ⁢                                                  ∂                          y                                                                                                      ⁢                                      |                    0                                    ⁢                                                            μ                      xy                                              (                        2                        )                                                              +                                                                  1                        2                                            ⁢                                                                                                    ∂                            2                                                    ⁢                          I                                                                          ∂                                                      y                            2                                                                                                                                ⁢                                      |                    0                                    ⁢                                                            μ                      yy                                              (                        2                        )                                                              +                    …                                                                                                            (          2          )                    Here μ(0), μ(1) and μ(2) are the zeroth, first and second moments of g, respectively:μ(0)=∫g(r,s)drds,  (3)μu(1)=∫ug(r,s)drds,  (4)μuv(2)=∫uvg(r,s)drds,  (5)where u and v correspond to r or s as appropriate.
By convention, the filter is normalized to unit integral and zero mean, soμ(0)=1,  (6)μ(1)=0,  (7)and the second moment becomes the covariance matrix of the smoothing filterVuv=∫(u−μu(1))(v−μv(1))g(r,s)drds=μuv(2).  (8)By plugging Equations 6-8 into Equation 2, it becomes clear that the leading term contributing to the bias introduced by the filter is due to the sum of the products of the components of the Hessian curvature matrix of the image and the covariance matrix of the filter.
                              J          -          I                =                                            1              2                        ⁢                                                            ∂                  2                                ⁢                I                                            ∂                                  x                  2                                                              ⁢                      |            0                    ⁢                                    V              xx                        +                                                            ∂                  2                                ⁢                I                                                              ∂                  x                                ⁢                                  ∂                  y                                                              ⁢                      |            0                    ⁢                                    V              xy                        +                                          1                2                            ⁢                                                                    ∂                    2                                    ⁢                  I                                                  ∂                                      y                    2                                                                                ⁢                      |            0                    ⁢                                    V              yy                        +            …                                              (        9        )            
The same arguments that lead to Equation 9 can be applied to the gradients and higher derivatives of the image. The leading term in the bias is due to a sum of products of the components of derivative tensors of the image with the covariance matrix of the filter. The filter is therefore also responsible for distorting the shape of the image.
The relationship between the denoised image J, the image I and the filter g is a convolution (Equation 2). It is therefore useful to take advantage of the convolution theorem to express the relationship as a simple product in Fourier space.{tilde over (J)}(k)=Ĩ(k){tilde over (g)}(k),  (10)where k is the vector wave number in Fourier space.
Noise is conveniently expressed in terms of its power spectrum in Fourier space. Denoising, therefore, results in a wave-number by wave-number attenuation of the amplitude of the noise by the noise-reduction factorNRF=|{tilde over (g)}(k)|  (11)For example, for a Gaussian filter{tilde over (g)}(k)=exp(−σ2k2/2),  (12)where k=|k| is the length of the wave-number vector and σ is the standard deviation of the filter, which is related to its full width at half maximum byFWHM=√{square root over (8 ln 2)}σ≈2.35σ  (13)
In general, filters transmit the lower wave numbers (longer wavelengths) and attenuate the high wave numbers, which replaces the image intensity in each pixel with a weighted average over the image intensities in its neighborhood. The details of the weighting scheme depend on the filter, but the qualitative behavior is to transmit the lower wave numbers and attenuate the high wave numbers. The transition from the low k to the high k regimes should also be smooth, in order to avoid oscillatory behavior due to the Gibbs phenomenon.
“Deblurring” is the process by which one attempts to undo any physical blur caused by the optical system or by the intervening medium between the source of the image and the detector. Blur also often takes the form of a convolution, in this case between the underlying true image and a point-response function, which is the image that a point source extends on the focal plane. In Fourier space, the blurred image is a product of the true image and the point-response function{tilde over (B)}(k)=Ĩ(k){tilde over (p)}(k),  (14)and the deblurred image can be written as a filtering with the inverse point-response functionĨ(k)={tilde over (B)}(k)/{tilde over (p)}(k).  (15)
Unfortunately, this process amplifies noise because {tilde over (p)}(k) is small at high wave numbers, and dividing by it magnifies the noise component of {tilde over (B)}(k). For example, for a Gaussian point-response function (PRF)1/{tilde over (p)}(k)=exp(+σ2k2/2)  (16)
which grows exponentially with wave number k.
The point-response function is determined by the physics of the optical system and the intervening medium. It is therefore not possible to modify the inverse point-response function. The only freedom available in image reconstruction is to modify the denoising filter.
The reconstruction method need not be an explicit Fourier technique, and it may also apply different degrees of denoising in different parts of the image, but the same principles apply. At the base of all image reconstruction methods is a selection of an explicit or implicit denoising scheme to either reduce the initial noise in the image and/or to suppress the noise amplification resulting from deblurring.
Tschumperlé, in his papers “Curvature-Preserving Regularization of Multivalued Images Using PDE's” (in European Conference on Computer Vision, Part II, Leonardis A, Bischof H, Pinz A (eds.), LNCS 3952:295-307, 2006) and “Fast Anisotropic Smoothing of Multivalued Images Using Curvature-Preserving PDE's” (Intl. J. Compu. Vision, 68:65-82, 2006) has recently proposed using partial differential equations to effect curvature-preserving regularization of images. This approach uses nonlinear, anisotropic smoothing using a vector or tensor field defined on the image, which preserves image curvature along curves defined by the vector or tensor fields. Tschumperlé uses the term “curvature” to mean a different quantity than the Hessian curvature matrix referred to herein. Also, this approach requires nonlinearity and anisotropy, while the curvature-preserving filters described herein can, and often will, be both linear and isotropic.
The more sophisticated reconstruction methods are, in general, iterative and therefore slower. The preference for iterative methods is due to the deficiencies of noniterative methods. Where noniterative methods suppress noise, they blur images; where they deblur, they amplify the noise in images. Puetter, Gosnell & Yahil (2005, Annu. Rev. Astron. Astrophys, 43:139-194), which is incorporated herein by reference, provides a comprehensive review of digital image reconstruction.
In the final analysis, digital image reconstruction is a tradeoff between competing effects, and the choice of the optimal reconstruction method depends on the goals of the user. It is therefore not possible to evaluate the adequacy of prior art without specifying the user's goals. For example, reconstruction of a still image makes fewer demands on image reconstruction processes compared to rapid processing of video in real time.